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Abstract 

We use 1-dimensional numerical simulations to study spherical collapse in the f{R) gravity 
models. We include the nonlinear coupling of the gravitational potential to the scalar field in 
the theory and use a relaxation scheme to follow the collapse. We find an unusual enhancement 
in density near the virial radius which may provide observable tests of gravity. We also use the 
estimated collapse time to calculate the critical overdensity 5 C used in calculating the mass function 
and bias of halos. We find that analytical approximations previously used in the literature do not 
capture the complexity of nonlinear spherical collapse. 
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I. INTRODUCTION 



The energy contents of the universe pose an interesting puzzle, in that general relativity 
(GR) plus the Standard Model of particle physics can only account for about 4% of the 
energy density inferred from observations. By introducing dark matter and dark energy, 
which account for the remaining 96% of the total energy budget of the universe, cosmologists 
have been able to account for a wide range of observations, from the overall expansion 
of the universe to the large scale structure of the early and late universe [lj. The dark 
matter/dark energy scenario assumes the validity of GR at galactic and cosmological scales 
and introduces exotic components of matter and energy to account for observations. Since 
GR has not been tested independently on these scales, a natural alternative is that GR 
itself needs to be modified on large scales. Two classes of modified gravity (MG) models 
are higher dimensional scenarios such as the DGP model, and modifications to the Einstein- 
Hilbert action known as f(R) models 2|, |3|, |4j. By design, successful MG models are difficult 
to distinguishable from viable DE models against observations of the expansion history of 



the universe. However, in general they predict a different growth of perturbations w 
be tested using observations of large-scale structure (LSS) (5], [f], B^, Q 
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Recently the nonlinear regime of structure formation in MG theories has been explored 
through simulations and analytical studies. Both the DGP and f{R) models have a mech- 
anism that restores the theory to GR on small scales. Recent studies of f(R) theories have 
focused on the effects of the chameleon field which alters the dynamics of mass clustering in 
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high density environments such as galaxy halos. A series of papers 
explored the consequences of this evolution through simulations and comparison to analyt- 
ical predictions. Similar efforts have been made to study large scale structure formation in 
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The evolution of isolated spherical overdensities in an expanding universe provides a 
useful approximation for structure formation in the universe. Although spherical collapse 
is just one idealized model for the formation of structure, it captures many crucial features 
of realistic mass distributions. This has been demonstrated by its successful applications in 
predicting the halo mass function, halo bias and merger history in the ACDM cosmology. 
Spherical collapse is sensitive to the nature of gravity, matter and energy. 
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Spherical collapse in general relativity with cold dark matter and smooth dark energy is 
unique in several aspects. For a spherical shell enclosing a fixed mass M, its collapse rate 
does not rely on the environment nor the internal density profile. Furthermore, an initial 
tophat spherical region remains a tophat during the collapse. Finally, after virialization, we 
have the simple relation 2K + W = between the total kinetic energy K and the potential 
energy W. All these interesting properties rely on the r~ 2 behavior of the gravitational 
force. Modifications to general relativity often destroy the characteristic properties described 
above: spherical collapse can become dependent on environment and internal structure, 
rendering an initial tophat density profile non-tophat and changing the conversion efficiency 
from potential energy to kinetic energy. If gravity indeed deviates from GR, we expect that 
at least some of these modifications would survive in realistic galaxies and galaxy clusters 
and serve as tests of modified gravity. 

Due to the complicated behavior of gravity in f(R) models spherical collapse no longer 
has analytical solutions. The purpose of this paper is to study spherical collapse through 
1-dimensional numerical simulations. Schmidt et al. [23] have used large-scale simulations, 



coupled with analytical approximations for spherical collapse, to predict the halo mass func- 
tions, linear bias and density profiles for the f(R) model of Hu & Sawicki [22j and compared 
them to the standard model of cosmology - the AC DM model. In addition analytical 
calculations have been done in the two limiting cases of the f(R) model: 

1. The strong field regime, where f(R) behaves like ACDM, but with a larger Newton's 
constant (by a factor of 4/3). 

2. The weak field regime, where there is no observable difference from the Standard 
Model. 
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The results from these bounding situations have been compared to simulations by 
the observed differences have been discussed. Since the strength of gravity in f(R) gravity 
lies inside these two limiting reasonable expectation is that the evolved observable 

quantities should also lie within the limiting cases. We will explore the validity of such an 
assumption by performing a direct simulation of a spherical collapse of an isotropic object. 
As chameleon f(R) theories exhibit highly non-linear behavior and there exist coupled fields, 
it is worth checking through explicit calculation the naive expectation based on limiting 
cases. 
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As described in 



24j | for simulations of f(R) gravity, the solution for the potential driving 



the dynamics of the evolution is coupled with the solution for the scalar field 22J. In 
our case we deal with isotropic objects and thus have a one-dimensional system. In §11 we 
present the radial equation for the fn field. In §111 we describe the simulation scheme for 
numerically solving the aforementioned equation. §IV describes the results, focusing on the 
distinct features of spherical collapse in f(R) gravity, while §V connects our results to the 
mass function. We conclude in §VI. 



II. RADIAL EQUATIONS FOR SPHERICAL COLLAPSE 



In general f(R) models are a modification of the Einstein-Hilbert action of the form: 

-R + f{R) 



S 



d X\J — g 



2k 1 



(1) 



where R is the curvature and k 2 = 8nG. The particular form chosen by 22j is: 



f(R) = -m 



, c 1 (R/m 2 ) n 
c 2 (R/m 2 ) n + 1 



with 



2 k2 Po 
m = — — 



(8315Mpcy 2 ('^ 



(2) 
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The properties of the model are well described by the auxiliary field /r = . ACDM 
expansion history with a cosmological constant Qa is obtained if we set: 



C2 



We choose to work with n — 1, which leaves one free parameter /rq ~ ~ n \ 



ci / 12 



9 



(4) 



-n-l 



to parametrize the model. Following 24Jj which set up the 3D simulation framework for the 
f{R) chameleon model, we start with the trace of the Einstein equations in the quasi-static 
limit, and the Poisson equation: 

1 

: 3c" 2 
167rGp 



V 2 /^ = [SR(f R ) ~ 8nG5p] 



(5) 



>Sp-±6R(f R ) 



(6) 



For the purposes of numerical calculations we need to define relevant dimensionless quan- 
tities, and we switch to comoving coordinates. Thus we adopt the definition of code units 



f=— , t = tH , p = a 3 — , 
r a p 



(7) 
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where 



Po — Pc,0^M,0, Rq 



r H 
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p = a- 



(r H ) 



v = r H 



(8) 



and r is an appropriate length scale (used for example to define the size of the overdensity) . 
Bare symbols X are physical coordinates/quantities while symbols with tilde X are code 
quantities, symbols with bars X are average physical and symbols with both bar and tilde 
X are average code quantities. 



Equations ([5] and [6]) then become the following in code units (Eqns. 25,27 in 241]): 

5R 



ac 2 



M,0 



5R 
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where 



5=^ = 1 



(9) 
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The next step is to express 5R in terms of fa in code coordinates. We start with (Eqn 9 in 
3): 



R = SixGp 



M 
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This is the average curvature in the f(R) model. Thus: 

R _ gAM / J_ ^ ^A,0 
Ro PO \« 3 ^M,0 



We arrive at: 



From that we also have: 



Ro \a 3 ^m,o 



(12) 
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We also need the relation between /r and R. Using (Eqn. 12 in {24]), defining fn(a = 1) 
fno, and working in the case n = 1 we see that: 

.//,> fR(" 1 ^ 
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This leads to: 



Ro R(a = 1 
For further use we will also need: 
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We can now obtain for the perturbation in the Ricci scalar: 

5R = R-R = ^-(R -R) = ^-5R = 
Ro Ro 



= 3a 6 1 + 4 
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From now on we will use only code quantities and drop the tilde notation. Noting that 



let's consider the following substitution: 



Jr — 



he" 



(21) 



Expanding the Laplacian (in code units) we arrive at the following equation for /#: 
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Additionally we will convert it to a system of 2 first order ODEs using an auxiliary function: 



d d 

y = 7r eU = eU ir u = eV 

Or or 



So the system with explicit r dependence looks like: 



(23) 



u\r) 
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r fl 



M,0 



o 3 / 1 + 4 ^ e _ u(r ) /2 _ _ 5(r) 



(24) 
(25) 



A detailed description of the relaxation method used to solve this system is presented in 
Appendix A. 

Considering that we can approximate 7/j ~ r^, where rjj is the upper boundary of 
integration in code coordinates, we can impose the boundary condition that /i?(r/j) = Jr 
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which translates to u(rfi) = ln(r/j). As the relaxation scheme requires 2 boundary points we 
will impose a condition on the inner boundary. We do not know the solution in the center 
of a collapsing isotropic mass distribution. What we know is that it has to be symmetrical. 
We also expect the solution in the center to be screened from the solution outside of the 
sphere by the thin screen (that is we expect to have chameleon effect). This means that 
very close to the center we expect to have behavior very similar to that of a homogeneous 
universe with that average density - but that would be a constant solution and thus zero 
derivative f' R (0) = 0. 



III. SIMULATION SCHEME 



To obtain the time evolution in the simulation, at each time step we proceed as follows: 

1. Given an initial density profile (from the previous time step) we compute the cor- 
responding solution for the fn field. Under the quasi-static approximation that we 
adopt, the gravity field is completely determined by the density distribution at the 
same epoch. 

2. This allows us to compute the solution for the Newtonian potential that drives the 
dynamics. 

3. The mass shells are then moved according to the dynamics equations 



da aa l 



and 

dp V0 



(27) 



da a ' 

where a = a" 1 / 2 jOj^o + Vl\fia 3 , as we tune the expansion of the universe to be the 
same as in AC DM. 

4. After the particles (shells) have been moved we can compute the new density distri- 
bution and proceed to a new time step thus closing the cycle. 
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FIG. 1: (Color online) Fractional residuals in the solution for fn at the final step of the relaxation process 
as a function of the time step. 

A. Code Tests 

An important issue in the case of numerical simulations is testing the code for stability 
and accuracy. The following have been checked: 



1. Self consistency: as per 24] we can start with an analytical function for fn(r). This 
can be analytically solved to obtain a corresponding density distribution. Now we can 
plug that density distribution in the numerical code and check how well the obtained 
solution reproduces the original analytical function. We observe deviations of the order 
less than 10~ 7 . 

2. During the relaxation scheme a measure of our accuracy is the residual relative size 
of the elements in the vector b (Eq j50|) as compared to the size of the corresponding 

elements of the solution — This tells us how much we need to correct the solution 

yfc 

obtained by the previous step in the relaxation. A sample plot of these as a function 
of time iteration step is provided at FigJTJ 

3. We have also checked the simulation for analytically solvable dynamics. In particular 



we observe that for simulations of 20000 time iterations we recover the analytical 
results (turn around radius, virial radius, collapse time) for the Einstein-De Sitter 
(matter dominated Q m = 1) model within 0.3%. This is the dominant numerical 
error. 



IV. VIRILIZATION IN f(R) GRAVITY 

An important question we need to address is how to identify the epoch at which a 
collapsing object reaches its virial radius. In the cases of Einstein-De Sitter (O^o = 1) and 
AC DM universes we have analytic solutions [23J] (Appendix A), but that is not so in the 
case of f(R) gravity. What we do employ is a step by step calculation of the Virial condition. 
Let's look at energy conservation. The total velocity v l = dx/dt = d(ar)/dt = dr + v, where 
x = ar is the physical distance, r is the comoving distance and v = ar is the proper peculiar 
velocity. The acceleration equation is 



d(av p ) d<fi 
dt dr 

On the other hand, v l satisfies another equation 



(28) 



u* = -^ ; 0* = ( f ) --adr 2 (29) 
dx ' y y 2 v ; 

Multiplying v l to both sides and integrating over t, we obtain the familiar energy conservation 

i(w*) 2 + 0* = constant (30) 

For this reason, v l and <p l are the relevant quantities for the virial theorem. Multiply Eq. 
[29] by x, we obtain 

_ M . _ -J*. (31) 

This equation is satisfied at all time. After virialization, we then take the average of the 
above equation for all particles. Now the velocity of particles is random (no correlation with 
x), so we have (xv 1 ) = 0. Then 

((^) 2 ) = <*^> (32) 

An equivalent expression, which can be applied straightforwardly, is 

r a At 

.t\2. 



2K = / (vydM = I x^dM . (33) 



Here, K is the total kinetic energy. The integral is over the region of mass M. This is the 
general expression of the virial theorem. One can check in the case of Newtonian gravity, 
0* oc M/x and J xdft/dxdM = — J § l dM = —W, this reduces to familiar form of the 
virial theorem, 2K + W = 0. See discussion in (30] for the pitfalls of using the Newtonian 
potential energy rather than the RHS of Eqn. 29 for modified gravity or general quintessence 
models. For the purposes of the simulation we need to express the above formulae in terms 
of comoving code coordinates given by Eqs. [7] and |HJ It is straightforwards to obtain: 

K = rl [ dM(ra+ 1 ^) (34) 



a 



for the kinetic term, and: 



W rr, / dMr { H$-^-raa\ (35) 

for the potential term. These are subsequently discretized and the sum is over the region 
with relevant mass. 

As expected in the case of GR (EDS and ACDM) the epoch at which the sum of these 
wo terms is zero coincides with the analytically predicted epoch of reaching the virial radius 



2J,|36|: 



Peff 2Q A (3g) 



' (1 + F)p m (l + F)n m a~ 3 (l + 5) 
2s - 1 

where s = t v /tta is the ratio of the virial radius and the maximal radius at turn-around. 
All relevant quantities are defined at turn-around. 

In our simulations we observe that the difference between the analytical result and our 
evaluation is of order 10 -5 for 20,000 time steps. We expect a similar level of accuracy to 
hold in the case of f(R) modifications. Thus we define the epoch of achieving virial radius in 
f(R) by the moment when this sum becomes zero during simulations. As per Fig. |2] observe 
that there are two moments when this condition is satisfied. We are obviously interested in 
the second one, which occurs after passing the turnaround point. 
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FIG. 2: (Color online) Evolution of the virial term. Observe that there are two points where it crosses zero. 
We are interested in the one that happens after turnaround. 

A. Density Enhancement at the Virial Radius 

For the study of spherical collapse in ordinary GR an initial top-hat density distributions 
is very convenient as it remains a top-hat during collapse. (In other words a top-hat is 
a Green's function for the spherical collapse evolution operator in GR). This allows for a 
straightforward definition of the key variables for spherical collapse - in particular 5 C and 
A vir . This is not the case for modified gravity theories. Unfortunately we do not know what 
profile would be the analogous Green's function. We can still compute the evolution of an 
initial top-hat distribution and try to compute 5 C and A vir in a similar way. 

We find that at the outer edge of the initial distribution the density becomes very large. 
This can lead to observable signatures for cluster halos, e.g. in weak lensing mass profiles. 
This effect appears qualitatively consistent with arising from chameleon screening. As the 
universe expands the size of the background fa field increases and we approach the high field 
limit of the f(R) theory - where it behaves as GR with enhanced Newton's constant. This 
means that the outer edge collapses faster than it would do in regular GR. The inside of the 
collapsing object, though, is under the effect of the chameleon and the solution for the 
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field becomes much smaller and thus the collapse slows down to approach the one in GR with 
Newton's constant. This makes the edge more and more dense as compared to the inside of 
the object, and this creates a positive feedback. The higher the edge density the stronger 
its screening effect and the inside slows down even further thus enhancing the accumulation 
of matter at the edge. However it has been suggested (Fabian Schmidt, private discussion) 
that the density enhancement we observe occurs for f(R) models with parameters that do 
not involve chameleon screening. A detailed study of environment and the nonlinearity of 
the theory is required to address the origin of the density enhancement. 

There is an interesting possibility (Justin Khoury, private communication): this density 
enhancement at the edge can separate the inside and the outside with an underdensity. We 
actually observe that the solution for starts exhibiting that kind of behavior in the very 
late stages of the collapse, but it is very close to the epoch of reaching virial radius so the 
effect is not observable in the density profiles. Clearly there are several issues in the late 
stages of spherical collapse that merit further study. 
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FIG. 3: (Color online) Comparison of the density profiles at the epoch of virialization for f(R) gravity (blue) 
and GR (red). In each case the starting profile (Gaussian smoothed tophat) and mass (1.5 x 10 14 M Q ) are 
the same. Virialization is reached at different epochs. 
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Comoving code coordinate r 

FIG. 4: (Color online) Smoothed density profiles with Gaussians with different dispersion. 

The edge effect leads to difficulties in the numerical integration for the initial top-hat. We 
therefore smooth out the edge with a Gaussian, which reduces the strength of the positive 
feedback and allows for stable evolution of the collapse. This smoothing is applied only to 
the original profile (step 1 of the simulation). The complicated part of that approach is that, 
as Birkhoff's theorem is not satisfied in f(R) models of gravity, the end result significantly 
depends on the environment, and in particular, what smoothing is used. 

The smoothed profile has one parameter, the dispersion of the Gaussian, which allows 
for controlling how close we are to a pure top-hat density distribution. The profile is: 

S(r) = MHe(r) - He(r - r TH )) + 5 in He(r - r TH )e^- r ^ (37) 

where He is the Heavyside step function. 

Even with smoothing, though, the code becomes unstable beyond reaching the epoch of 
the virial radius, preventing us from achieving collapse to singularity - the epoch of collapse. 
In Fig. [3] we show a comparison of the density profiles at virialization between f(R) and 
ACDM. In each case the starting profile and mass (1.5 x 1O 14 M ) are the same. They 
achieve virialization at different epochs. 
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FIG. 5: (Color online) Velocity ratio computed shell by shell and normalized by the physical position of 
the shells. Shell number 100 represents the edge of the top-hat part of the initial overdensity. 

V. IMPLICATIONS FOR THE MASS FUNCTION: THE COLLAPSE THRESH- 
OLD 5 r 



23j dealt with spherical collapse in an analytical way (Appendix A) by solving the 2 
limiting cases for the strength of the effective Newton's constant in the f(R) model of 
gravity. The prediction in the end is that the fundamental quantities should lie inside of 
the region bound by the values of the 2 limiting cases. In particular they are identified by 
the value of the parameter F which governs the strength of the effective Newton's constant. 
Regular GR corresponds to F = and the strong field limit to F — 1/3. For Qm,o — 0.24 
these imply 5 C = 1.673 for F = 0, and 5 C = 1.692 for F — 1/3. 

Computing S c is straightforward in ACDM with GR. For a given starting epoch ai n we 
need to find an initial overdensity S^gr-, which would collapse to a singularity at the present 
time. Then we just need to evolve that initial overdensity to present time via the linear 
growth factor. In the case of ACDM this can be performed analytically ({3?]] Appendix A). 
In the case of f(R) modified gravity there are 2 complications. 

1. The linear growth factor is scale dependent. 
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2. Our simulation allows us to only reach the epoch of reaching the virial radius and not 
the epoch of collapse. 

Resolving the first issue is not a complicated task. The solution is to go to Fourier space 
and convolve the linear growth factor at the epoch of collapse (normalized with the growth 
factor at the initial epoch) with the Fourier image of a top-hat function. After that, we need 
to Fourier transform back to physical space, which is greatly simplified, as we are interested 
only at the value at r = 0, and sums up to the evaluation of an integral. 
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FIG. 6: (Color online) S c as a function of field strength /rq for mass (1.5 x 10 14 M Q ). The results are color 
coded to represent the dispersion of the smoothing Gaussian. 

Next we deal with the problem of estimating the collapse epoch. The numerical issues we 
have with the development of a density spike at the edge require the use of approximations. 
First we study the effect of the environment due to the invalidity of Birkhoff's theorem. 
Recall that if Birkhoff's theorem is valid in a gravitational theory then the behavior of a 
shell depends only on the mass inside the ball enclosed by that shell. In our case that is 
not correct and shells are influenced by what is outside - the environment. We utilize a 
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sequence of initial profiles, each of which has a pure top-hat part and then is smoothed with 
a Gaussian with varying dispersions FigJH that bring us closer and closer to a pure top-hat 
distribution. This way we can study the trend of changes in the top-hat part of the initial 
profile when approaching a pure top-hat overdensity. 

In order to estimate the collapse epoch, we first note that the time/scale factor between 
the epoch of reaching the virial radius and the epoch of collapse to singularity is a small 
fraction of the total time/scale factor in the evolution of the spherical object. Thus we 
will assume that if an object in f(R) gravity achieves its virial radius at the same epoch 
as a corresponding object in ACDM does, then these two objects should reach collapse to 
a singularity at approximately the same epoch as well. Then we can make an estimate of 
how wrong we are in this prediction. The task of finding 5 C then is moved to finding the 
initial overdensity in f(R), which would reach its virial radius at the same epoch at which 
a corresponding object in GR does. In addition we require that the GR object collapses to 
a singularity at the present epoch. 

What is left is estimating the error of this calculation. One way to approach the problem 
is to look at the radial velocity field of the evolving shells in our simulation and compare 
them between the corresponding objects in f(R) and ACDM. In Fig. [5] we show the 
velocity ratio computed shell by shell and normalized by the physical position of the shells 
(the corresponding f(R) and ACDM have different mass and different size at the epoch 
of achieving the virial radius). The different colors correspond to the different smoothing 
factors we have introduced as way to approach a pure top-hat distribution. In our simulation 
we have chosen shell number 100 to represent the edge of the top-hat part of the initial 
overdensity. As we can see the normalized velocity ratio remains within 5% of unity at the 
edge of the top-hat, which suggests that a good estimate of our error would be of the same 
order. 

Another way to approach the issue is to vary the initial overdensity and look at how much 
it changes the epoch of achieving the virial radius and compare with the expected epoch 
of collapse. In particular, values for the initial overdensity, that have epoch of virial radius 
close or beyond the expected epoch of collapse, set a bound on our error. We found that 
this also puts a hard error bar of 5%, which is what we finally used in our calculation. 
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VI. RESULTS AND DISCUSSION 



The main results of this work is presented in Figs. [3]and|6j The development of an excess 
overdensity at the edge of spherical halos in f(R) gravity is shown in Figs. [3J While we 
have not carried out detailed studies of realistic mass profiles, the results suggest that the 
region around the virial radii of cluster halos may contain signatures of f(R)— type theories 
of gravity. 

For the collapse threshold 5 C , our results are shown in Fig. |6j We tested the following 
conjectures: 

1. In the weak field regime our calculations should approach the result for regular strength 
AC DM (F = 0). 



2. In the strong field regime the result must approach the values predicted in [23[ for 
F — 1/3. This behavior is not guaranteed. We know, for example, that in this limit 
Eqfl2l is not valid. 

We find a significant dependence of the values of 5 C on the field strength, and particularly 
in the physically interesting region around fno = — 10~ 6 - currently close to the upper bound 
permissible by observations or theoretical considerations. We observe a strong environmen- 
tal dependence with a significant trend: when reducing the dispersion of the smoothing 
Gaussian (and thus approaching pure top-hat distribution) we deviate further away from 
the analytical prediction in |23J . This result shows that the non-linear chameleon properties 
of the f(R) models strongly affect its behavior; thus analytical approximations based on 
linear predictions should be viewed as simple guidelines. However our error bars are still 
large, a more careful study is needed to make any definitive statements about S c . 

Studying smoothed top-hat initial profiles and obtaining estimates for S c is the first step 
in studying spherical halos in f(R) gravity. Further work is needed in understanding realistic 
halos with differing masses and environment. It would also be interesting in future work to 
study the abundance and clustering properties of halos: mass function and halo bias. 
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Martino, and Ravi Sheth for helpful discussions. We are grateful to Fabian Schmidt for 
comments on the manuscript and to him and Lucas Lombriser for sharing the results of 
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APPENDIX A: RELAXATION SCHEME FOR SOLVING THE SYSTEM OF 
NONLINEAR ODES 



As discussed in [22| (p. 8) the primary equation we need to solve (EqEland subsequently 
EqJiMl and Eq j2"o"j) is non-linear and cannot be solved as an initial value problem as the 
homogeneous equation has exponentially growing and decaying Yukawa solutions e^ ±ar '/r. 
Initial- value integrators have numerical errors that would stimulate the positive exponential, 
whereas relaxation methods avoid this problem by enforcing the outer boundary at every 
step. So we also employ a relaxation method for solving 2-point boundary problems in ODE. 



We employ a Newton's method 35( with dynamical allocation of the mesh grid. The mesh 
allocation function is taken to be logarithmic with its higher density at the origin, which is 
the primary region of interest and where we expect the solution to be more rapidly changing. 
As a guess solution for each step we utilize the relaxed solution of the previous step, while 
for the initial guess at the beginning of the simulation we use a linear solution. Generally if 
we have a system of discretized first order ordinary differential equations in the form: 

= E fc = y k - y fc _i - (x k - x k ^)g k (x k , x k -i, Yk, Yk-i) k = 2..M, (38) 

where the index k spans the number of grid points 2..M, the vector E consists of the system 
of N discretized 1st order ODEs at each point (and has a total of N * M components - 
N * (M — 1) from differential equations and N from boundary conditions). Ei and Em+i 
describe the boundary conditions. So a Taylor expansion with respect to small changes Av/c 
looks like: 

E k {y k + Ay k , y fc _x + Ay fc _i) « (39) 

, . x -* dE k ^ — \ <9Efc 
w E k (y k , y fc _i) + } &y n ,k-i + > v -k &yn,k, 

For a solution we want the updated value E fc (y fc + Ay^y^ + Ay fc _x) to be zero, which 
sets up a matrix equation: 

TV 2N 

^ S jtr Ayn,k-l + ^2 S 0,n^yn,k = ~Ej,k, (40) 
n=l n=N+l 

where 

_ dE j)k _ dE jyk 

oy n ,k-i oy n ,k 
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and the quantity Sj >n is a N x 2N matrix at each point. Analogously we obtain similar 
algebraic equations on the boundaries. Considering our problem we look at Eqs. (1241 and 
to obtain (after discretization and using the appropriate variables): 

r k + r fc _i \ 1 n Mfi 



Eih = (y k - Vk-i) ~ (r k - r fc _ 



fa «c : 



(42) 



a s — + 4 



^A,0 



lr k + r fc _i _ ( '- fc +r fc _ 1 \ r +r 



E- 



2,k 



2 



[u k - M fc _iJ - (r fc - r fc _ij | I < 

<Sl,l,fc = — 1) $l,3,k — 1 



Sl,2,k — SiA k 



4 



(43) 
(44) 

(45) 



1 



+ 4 



Af,0 



^A,0 ^ ^ / rfc + rfc - 1 e ~( rfc+ 4 fc ~ 1 ) 



5- 



r k - r fc _x \ _ f ^vi) 

1 v 2 . 



5. 



2,2,fc 



2,l,fc — 52,3,A; 



r k - Tk-A f y k + y k -i\ _ ( -k-r k -i ) 



5- 



2,4,Jfc 



+ 1 



2 y V 2 

The boundary conditions are also easily translated in terms of the relaxation scheme. 



E- 



2,0 



(46) 

(47) 
(48) 

(49) 



The notation here is probably a bit confusing. The quantity E is a vector which consists 
consecutively of 2 elements per M — 1 grid points. In addition there is one element each 
at the beginning and the end that correspond to the boundary conditions. Thus the total 
length of E is NM, while the matrix S has NM x NM elements. The task of relaxing the 
solution at each step of the scheme requires solving the matrix equation: 



Sb = E 



(50) 



The vector b contains the updates Ay^. For a grid of 1000 points our equation requires a 
matrix of derivatives of size 2000x2000 elements. Fortunately it is sparsely populated and 
as such can be represented by a sparse array structure. This allows for the use of methods 
particularly designed for solving such systems, like Krylov's method, which we employ. 
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APPENDIX B: SPHERICAL COLLAPSE IN ACDM 



The evolution and collapse of spherical overdensities have been useful for modeling the 
formation of galaxy and cluster halos. In GR the problem can be approached analytically 
and we will outline the derivation presented e.g. in [23]. We start with the nonlinear and 
Euler equation for a non-relativistic pressureless fluid in comoving coordinates: 

^ + -V-(l + 5)v = 

at a 

dv 1 1 

- + -(v-V)v + ffv = V0, (51) 

at a a 

where a(t) is the expansion scale factor, H(t) = a/ a, and is the "Newtonian" potential. 
These equations continue to be valid for modifications of gravity that remain a metric theory 



381 ] . These can now be joined together to form a second order equation for S. 



d 2 5 | 2H d5 1 d 2 (l + S)vhP _ V ■ (1 + £)V0 
dt 2 dt a 2 dx l dx^ a? 

Solving this equation requires information about the velocity and potential fields. In the 
case of a spherical top-hat distribution, to preserve the top-hat distribution, the velocity 
field must take the form v = A(t)r to have a spatially constant divergence. Its amplitude is 
related to the top-hat density perturbation through the continuity equation: 

S+-(1 + 5)A = (53) 

(X 

This leads to: 

dV^' 10 . 2 4 2 5 2 

„ . = 12A 2 = -a 2 - — . (54) 

d&dx? 3 (1 + 5) 2 v ; 

Substituting the above relation, the equation for the evolution of S becomes: 

PS 86 4 5 2 (1 + ^ 

w + 2H d-t-3jrTT) = ^^ v ^ (55) 

which is completed by the Poisson's equation for the potential: 

V 2 = 47rGa 2 Sp m . (56) 

It is common to express spherical collapse through the evolution of the radius of the top-hat. 
For that we use mass conservation: 

47T 

M = — r 3 p m (l + 6) = const. (57) 
o 
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to obtain the following relation: 



- = H 2 + H--^ (58) 
r 3a 2 



Expressing derivatives in terms of scale factor ' = d/dlna, with the useful substitution: 



w = , (59) 



and using Poisson's equation we obtain: 



, H 1 n m a - 2tt A 1 tt m a a 
w + -T7 W — 5 7^ — w ~~ TTT^ 5 T^TK h WjA, (60) 



where 

3 



A= — + (61) 

\1 + waif a J 

In these coordinates collapse occurs when w = — — . The task of computing 5 C now reduces 
to the following: for a given a, find an initial overdensity 5{ such that the collapse occurs at 



a = 1. Then using the linear growth factor in AC DM (see for example 38]) we extrapolate 
Si to the present epoch to obtain 5 C as: 

5 c (r = 0) = A [ D ^ , > a= - L ^ (fc,a in )e ffcr | r= odfc (62) 
J D(k,a in ) 
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